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A procedural code (C++) simulation was developed to investigate potentials for 
mixture ratio control of pressure-fed spacecraft rocket propulsion systems by 
measuring propellant flows, tank liquid quantities, or both, and using feedback 
from these measurements to adjust propellant tank pressures to set the correct 
operating mixture ratio for minimum propellant residuals. The pressurization 
system eliminated mechanical regulators in favor of a computer-controlled, servo- 
driven throttling valve. We found that a quasi-steady-state simulation (pressure 
and flow transients in the pressurization systems resulting from changes in flow 
control valve position are ignored) is adequate for this purpose. Monte-Carlo 
methods are used to obtain simulated statistics on propellant depletion. Mixture 
ratio control algorithms based on proportional-integral-differential (PID) controller 
methods were developed. These algorithms actually set target tank pressures; the 
tank pressures are controlled by another PID controller. Simulation indicates this 
approach can provide reductions in residual propellants. 


I. Introduction and Purpose 

One of the least reliable components in a pressure-fed rocket propulsion system is the mechanical 
regulator that controls tank pressurization. In addition to reliability issues, regulators are difficult to set, 
and uncertainties are associated with the delivered pressure. Current spacecraft practice uses separate 
regulators for the oxidizer and fuel to prevent accidental mixing of propellants in the pressurization gas 
system during long planetary coast periods. This introduces an additional problem: differences in the 
output pressures of the two regulators can aggravate operating mixture ratio differences and increase the 
residual propellant component of inert mass. Allowance for residual propellant contribution to inert mass 
must be enough that the probability of depletion of propellant prior to mission completion is very small. 

The reported investigation was undertaken to determine the performance potential for active 
mixture ratio control systems for representative planetary spacecraft, thus to assess whether this 
technology is a promising area for technology advancement investment. Most if not all current spacecraft 
propulsion systems are pressure-fed, and this investigation was limited to pressure-fed systems. 

Most current pump-fed engines are capable of changing mixture ratio over a modest range, and 
some launch vehicles use this capability for mixture ratio control. The main need for launch vehicle 
applications is better tank liquid quantity and flow measurement instruments. 
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A note on units: We have taken the liberty of using SI units for all parameters in this paper 
except pressures. Most propulsion analysts are familiar with system pressure characteristics quoted in 
psia, so we follow that practice; kPa values are given in parentheses in some cases. Simulation 
calculations all used SI units. 


II. Initial Simulations 

Simulations of tank pressure control began a few years ago at Space America, during tests of a 
pressure-fed LOX-kerosene engine. Test stand experience with propellant pressurization regulators was 
less than satisfactory. An initial system using on-off solenoid valves was targeted. The method was 
simulated prior to implementation, to determine practicality and develop controller algorithms. A 
proportional-differential controller was used, with input from tank pressure transducers and the output 
digitized to an on/off time ratio. 

The test stand computer control cycle was 50 milliseconds. The valves responded adequately, but 
the minimum pressurant flow bit (i.e. the minimum amount of gas delivered by one minimum-duration 
valve on-off cycle), with helium bottles at full pressure and minimum ullage in the tanks, led to pressure 
variance in ullage pressure about 5 psi. However, once propellant flow started, regulation was more 
accurate and repeatable than mechanical regulators. 

This experience led to a concept using multi-tum throttling valves driven by “go-to” digital servo 
steppers. It was simulated, but never implemented due to lack of funds. That simulation was the starting 
point for the present effort. 


III. Mixture Ratio Control Approaches 

When separate active controllers are used for the two propellant tanks, the opportunity exists to 
adjust tank pressures to control the system operating mixture ratio. Variance in operating mixture ratio 
due to uncertainties in pressures and pressure drops are only a few percent. Tank pressure adjustments to 
control these are also only a few percent. If these variances are not controlled, residual propellant due to 
one propellant depleting sooner than the other is also only a few percent of loaded propellant, but in a 
typical propulsion system, that represents a significant fraction of propulsion system inert mass. 

In order to effect propellant simultaneous depletion by controlling mixture ratio, accurate 
measurements are necessary. We considered three operating scenarios: 

1. Amounts of propellants loaded are precisely known, for example, because the propellants are 
storable and are loaded in a propellant servicing facility where the spacecraft is weighed during loading. 
Flow measurement devices (flowmeters) are available for measuring propellant flows during system 
operation. 

2. Amounts of propellants loaded and/or available are not precisely known, for example, because 
they are loaded on the pad. Also, one or both propellants may be cryogenic and some venting of boiloff 
may have occurred during a coast period prior to propulsion system use. Accurate measurement of liquid 
quantity in each tank is necessary. Liquid quantity measurements, either continuous or periodic (these 
cases are different!), are to be used for mixture ratio control. 

3. As in case 2, propellants available are not accurately known. Liquid quantity measurements 
and flow measurements are to be used. 


IV. Simulation Target 

The system simulated was a conventional spacecraft propulsion system. The schematic used for 
the simulation is shown in Figure 1. The figure shows half of the propulsion system, i.e. the oxidizer 
side. The fuel side is identical. The heat exchanger is optional and can be left out of the simulation. The 
purpose of simulating a heat exchanger was to evaluate the benefits of storing the helium pressurant cold 
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(to reduce volume and mass of the helium tanks) and warming it before it enters the propellant tank. 
Cavitating venturis are also an optional feature. These are used to make a system insensitive to pressure 
drops downstream of the cavitating venturis. These devices when operating properly, i.e. cavitating, 
exhibit mass flow that is a function only of the upstream pressure. Downstream pressures adjust to the 
flow set by the cavitating venturis. 

Typical systems operate with tank pressures 250 to 300 psia (1720 to 2070 kPa), combustion 
chamber pressures about 140 psia (965 kPa), thruster injector pressure drop about 1 10 psia (760 kPa), and 
cavitating venturi pressure drop (if used) about 35 psi (240 kPa). Other propellant flow pressure drops are 
small. 

Helium tank initial pressures are typically 4500 psia (31 MPa), and at the end of propellant 
expulsion about 450 psia (3.1 MPa). 
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Figure 1: System Schematic Showing Oxidizer Side (Fuel Side is the Same) 
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The simulation considered the following pressurization system segments: (1) Source tank, (2) 
Line 1, from the tank to the control valve, (3) control valve, (4) Line 2, from the control valve to the heat 
exchanger, (5) Heat exchanger, (6) Line 3, from heat exchanger to propellant tank, (7) Propellant tank 
ullage, which varies in volume as propellant is expelled, (8) propellant tank heat capacity and heat 
exchange with the pressurant gas in the ullage volume (heat exchange with the propellant is also 
estimated, but propellant evaporation into the ullage volume was not simulated; this has essentially no 
effect on results for storable propellants that were the subject of this study), (9) liquid feed system 
pressure drop, (10) injector pressure drop, and (11) thruster chamber pressure. The latter was modeled as 
proportional to total propellant flow. 

The controller simulated was assumed to be a throttling valve such as a needle valve, with 
roughly linear characteristics over the control range of travel. The valve is driven by a digital servo with 
rotary motion. We assumed valve travel of three complete rotations and a servo capable of 1 -degree 
steps. This gave 1080 steps total, which we rounded to 1000. The control algorithms do not depend on 
any particular valve calibration. 
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V. Simulation Methods 


Dynamic Simulation - We refer to dynamic simulation as a method that attempts to keep track of 
changing pressures, pressure drops, temperatures, and mass flows in and out as a function of time and 
controller settings. Mass flows in and out of each segment of the system are in general not equal. The 
first attempt at dynamic simulation tried to evaluate flow from the source tank through the system. Flow 
at the control valve was calculated first (it’s usually choked), and calculations were then propagated 
sequentially upstream and downstream. This method was not satisfactory as it exhibited calculation 
instabilities and did not match conditions at the interfaces between system segments. Typical results are 
shown in Figure 2. The left-hand plots are pressures, the upper right is temperatures, and the lower right 
is flows. The upper left shows instabilities with a large time step. At the lower right, the instabilities are 
reduced by a much smaller time step. The stair-step shape of the curve results from stepwise opening of 
the control valve. The lower right shows that flows do not converge to a common value as they should. 
Interface conditions are not matched. 
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The second approach used a simultaneous solution based on first-order time derivatives and a 
Taylor’s expansion. In this way, all the interactions among variables were solved jointly. If one 
considers, for example, the outflow from the storage tank: it is determined by the (small) pressure drop 
from the storage tank into Line 1: 

m-Cp A/sqrt(RT) 
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where C is the flow coefficient, and p and T are upstream values. R is the gas constant, and A is flow 
area. The variables are C, p, and T. The Taylor’s expansion is 

dm = (Dm/DC) SC + (Um/Up) dp + (Um/UT) ST 

Analogous expressions are written for all of the variables, and the partial derivatives are calculated. 
(They may be calculated from analytical expressions or derived in a computer program from the 
expressions for the variables by numerical differentiation; we chose to use analytical expressions.) The 
Taylor’s expansion expressions are then placed in a linear algebraic matrix as illustrated in Figure 3. The 
matrix is inverted to obtain a simultaneous solution for all the changes (5’s) in parameters, for a small 
time step. The process is repeated as a numerical integration. The test matrix shown in Figure 3 
simulated the helium storage tank, Line 1, and the flow control valve. (Items in red were colored to aid in 
correcting errors.) “T” items are Taylor’s expansions too lengthy to write into the matrix diagram 
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Figure 3: Test Matrix for Partial Simulation 


Typical test results are shown in Figure 4, and showed (a) that the method works, i.e. initial 
values of parameters change to a steady-state in which flows and pressures are in balance; mass changes 
in Line 1 go to zero as steady state is reached; (b) the time to reach steady state is very short, typically 10 
microseconds. Time steps about 1 microsecond were needed. In the simulation, the cycle time between 
helium control valve steps was 50 ms (50,000 pis). Therefore, we decided that simulating transients in the 
helium pressurant delivery system was unnecessary as well as impractical. A quasi-steady state solution 
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method was adopted, in which it is assumed that helium flow is set by the control valve and is the same at 
any particular instant everywhere in the helium delivery system. 
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Figure 4: Typical Results of Dynamic Simulation 


Ouasi-Steadv-State Simulation - Flow is set by the control valve position. A pressure iteration is 
needed for the pressure drop through the control valve, which normally will be choked. The flow in Line 
1 is always subsonic. The flow in Line 2 may be partially supersonic and shock down to subsonic flow. 
Constraints are that the static pressure at Line 3 exit into the propellant tank must be equal to the current 
ullage pressure (assuming the flow here is subsonic). Pressure downstream of the valve in Line 2 is 
iterated until this constraint is satisfied. Simulation time steps were half the assumed control cycle time, 
i.e. 25 ms. Figure 5 shows an early test of this simulation, in which LOX tank ullage pressurization is 
simulated. The pressurization algorithm steps the valve open in discrete steps until the ullage pressure is 
about half the target value; then begins stepping the valve closed. When the ullage pressure reaches 90% 
of the target value, control authority is transferred to the run pressure control algorithm. 

VI. Pressure Control Algorithms 

Two pressure control algorithms were tried: a predictor-corrector and a proportional-differential 
(PD) controller. The predictor-corrector attempted to determine the valve position that would result in the 
desired ullage pressure as a quasi-steady-state, based on measured values of current pressures, valve 
position, and rates of change of pressures. It then drove the valve to that position. The PD controller 
calculated a value based on gains applied to the pressure error and the rate of change of pressure. This 
value was converted to a valve step increment (amount of rotational motion), which was then executed. 
The valve step allowed for a single control cycle was limited to a maximum value to avoid excessive 
valve motion. Both worked; the PD controller was simpler and a little more stable, with less valve 
dithering. 
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Figure 5: Early Test of Quasi-Steady-State Simulation 


VII. Mixture Ratio Control Measurement Methods 

Mixture ratio control feedback may be provided by propellant flowmeters, by tank quantity liquid 
measurement, or by both. The NASA Marshall Space Flight Center is performing R&D on new 
technology for these purposes. For an active control system to exceed the performance of simple 
calibration requires that the measurements exhibit good accuracy. Historically, flow and quantity 
measurement methods have not offered enough accuracy to give an expectation of significant 
performance improvement. 

A balanced flowmeter 1 is being developed to measure liquid mass and volume flow. This is a 
pressure-drop device and does not have a rotor in the flow that could break up or otherwise fail in a 
hazardous manner. Two approaches are being taken to propellant tank liquid quantity measurement 2 ’ 3 ; 
one is a variation on the PVT method and the other is an optical device. Measurements during thrusting 
operations, when vehicle acceleration settles the propellants, could also rely on point liquid sensors. As 
noted below, a number of sensors would be required. 

VIII. Mixture Ratio Control Algorithms 

1. Loaded propellant quantities accurately known: This is a classical proportional, integral, 
differential (PID) controller problem. The error signal is difference between the current run mixture ratio 
as determined by flowmeters, and the known loaded mixture ratio. A PID controller applies a gain to the 
error signal and adjusts the tank “set” pressures [“Set” pressure is the commanded pressure input to the 
tank pressure controller.) based on the gain. The integrator in the algorithm accumulates the error and 
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applies a second gain. This causes the controller to correct for the small error the proportional factor must 
have to generate a corrective signal as well as prior accumulated error. The differential part of the 
controller applies negative feedback based on rate of change of error and stabilizes the controller. 

2. Current propellant quantities not accurately known (as in cases of boiloff, or use of one 
propellant as a monopropellant for attitude control): There are two sub-cases, 2a and 2b. In 2a, tank 
liquid quantity measurement is continuous or nearly so, and errors in the measurement do not vary in a 
random manner; that is, errors are systematic. In this case, an error signal can be derived from the rate of 
change in tank quantities mixture ratio. Driving this rate of change to zero means the tank quantities 
mixture ratio stays constant and both tanks will deplete simultaneously. In 2b, tank liquid quantity 
measurements exhibit errors that fluctuate in a random manner. This causes the 2a algorithm to become 
unstable because the differential part of the controller sees these random variations as large rates of 
change in tank mixture ratio, and outputs unacceptable fluctuations in tank “set” pressures. 

In the 2b case it is necessary to introduce some statistical smoothing to obtain a satisfactory 
algorithm. One way to do this is illustrated in Figure 6. Consecutive measurements are collected in 
circular buffers. Once enough are collected, regression lines are computed for rate of depletion of each 
propellant. These are extrapolated to zero propellant, which is the estimated depletion time. The 
difference in depletion times is the error signal. The controller attempts to drive this error to zero. 
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Figure 6: Regression Method for Using Tank Liquid Quantity Measurements 


3. This one works the same as the first one, except that current measured tank liquid quantities 
are used instead of initial values. Flowmeter readings provide the current run mixture ratio, and the 
difference between tank and run mixture ratios is the error signal. The continuing measurements of tank 
mixture ratio provide the integration function, so an integrator is not needed in the controller. If tank 
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liquid quantity measurements exhibit random fluctuations, current measured tank liquid quantities should 
be estimated from regression lines rather than from single measurements. 

IX. Monte-Carlo Simulations 

Monte-Carlo simulation enables estimates of the statistics of the output of a system based on 
predicted random fluctuations of the inputs. The method is conceptually simple: a sequence of random 
numbers is generated, usually varying from 0 to 1 (with median 0.5). The probability of occurrence of a 
random number within any particular range, say, 0.01 to 0.02, is equal to the probability of occurrence 
within any other equal range, say, 0.51 to 0.52, in the interval 0,1. In other words, the output of the 
random number generator is a “flat” (uniform) distribution. By virtue of the central limit theorem, 
random outputs of complicated physical systems vary according to the normal (Gaussian) distribution. 
The variances of inputs need to be converted from a uniform distribution to a normal distribution; 
therefore the random numbers in the sequence must be converted to a normal distribution. In this study 
we used a curve fit. 

Figure 7 shows the use of the normal distribution integral curve to obtain the equivalent variance 
in standard deviations (sigmas) from the output of the random number generator. The mean of the 
uniform distribution is 0.5. As an example, the deviation from the mean of the vertical arrow is about 
0.25 which corresponds to 0.68 sigma 
from the mean. If the random number 
generator outputs a value 0.75, the 
variation in sigmas for the simulation 
input is 0.68, and the 1-sigma value is 
multiplied by 0.68. Similarly, if the 
random number output is 0.1, the 
simulation input is varied by -1.2 
sigma. 

Variances were applied to 
inputs such as helium and propellant 
line pressure drops, thruster injector 
pressure drop, measurements of tank 
pressure, flowmeter measurements, 
and heat transfer in those cases where a 
heat exchanger was used. Early tests 
of the complete simulation, during 
which controller algorithms were 
developed, did not apply random errors The random number obtajned from the generator js input to this 

to flowmeters and tank liquid quantity curve to get the number of sigmas it represents. The data base 

measurements. Later simulations error value was taken as 1 sigma, and that value was multiplied 

included measurement errors. by the number of sigmas to get the varied value. The generator 

During preparation of this is 9 uer 'ed for each value from the data base, so the generator is 

paper, a better method of generating 9 ueried almost 40 times for each simulation run. 

normally distributed random numbers Figure 7; Use ofCurve FU tQ 0buUn Number ofSi 
was found. This uses a Box-Muller 

transform. The method directly transforms a uniform distribution to a normal distribution as summarized 
in Figure 8. Uniformly distributed random numbers are generated in pairs, and from these a pair of 
normally distributed numbers is generated. As noted in the Figure, this only requires two lines of code 
(three counting the extra call of the uniform random number generator. A preliminary test is shown in the 
Figure. Statistical tests on the normal distribution were not performed. 
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Figure 8: Box-Muller Transform Method 


ranxl = rangen (&seed1, &seed2, &seed3); 
ranx2 = rangen (&seed1, &seed2, &seed3); 
ranzl = sqrt(-2*log(ranx1))*cos(2*pi*ranx2); 
ranz2 = sqrt(-2*log(ranx1))*sin(2*pi*ranx2); 

• Rangen is the uniform distribution 
random number generator 

• Ranxl and ranx2 are the uniformly 
distributed pair 

• Ranzl and ranz2 are the normally 
distributed pair 


X. Simulation Flow Chart and Representative Results 
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A flow chart of the 
complete simulation is 
shown in Figure 9. The data 
base contains the inputs to 
the simulation, and can be 
edited on-line before 
launching the simulation or 
off-line. The Monte-Carlo 
driver loops around the 
entire simulation. Outputs 
include a time plot of one of 
the runs (100 Monte Carlo 
runs were usually used) and 
a summary of results of the 
Monte Carlo runs. 

An early simulation 
result is shown in Figure 10, 
no mixture ratio control on 
the left and with mixture 
ratio control on the right. 

The figure represents one time through the simulation. Without mixture ratio control, the ullage pressure 
controller maintains constant ullage pressure. The run mixture ratio is constant. Because the run mixture 
ratio is higher than the initially loaded mixture ratio, the mixture ratio remaining in the tanks decreases 
until oxidizer is depleted at tank mixture ratio zero. This case was a poorly calibrated system that had 
significant residuals in the no-control case. The curves on the right are with a PD controller (no 
integrator), perfect knowledge of initially loaded propellant, and perfect flowmeters. The mixture ratio 
controller alters ullage pressures, attempting to control the run mixture ratio to the initially loaded mixture 
ratio. Because the controller does not include an integrator, an error remains and gradually becomes 
larger. Tank pressures are adjusted by the controller until the pressures reach pre-set limits, and are then 


Figure 9: Simulation Logic Flow 
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constant. Even without an integrator in the controller, the residual propellant is reduced by 75%. An 
integrator was also tested, and in this case of perfect knowledge of tank liquid quantities and 
measurements, residual propellant decreased to near zero. 

Type 1, tank quantities known and flows continuously measured, no integ- 
rator in controller. When integrator was added; residuals decreased further. 
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Figure 10: Early Simulation Run Results Summary 
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